Transient dynamics of the Anderson impurity model out of equilibrium 
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We discuss the transient effects in the Anderson impurity model that occur when two fermionic 
continua with finite bandwidths are instantaneously coupled to a central level. We present results 
for the analytically solvable noninteracting resonant level system first and then consistently extend 
them to the interacting case using the conventional perturbation theory and recently developed 
nonequilibrium Monte Carlo simulation schemes. The main goal is to gain an understanding of 
the full time-dependent nonlinear current-voltage characteristics and the population probability of 
the central level. We find that, contrary to the steady state, the transient dynamics of the system 
depends sensitively on the bandwidth of the electrode material. 

PACS numbers: 73.63.Kv, 73.63.-b, 5.10.Ln 
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I. INTRODUCTION 

The Anderson impurity model (AIM) has been intro- 
duced in the early 1960s to describe conduction elec- 
trons interacting with a magnetic atom and since then 
continues to attract the attention of condensed matter 
physicists* 1 - Despite some notable exceptions^ for quite 
a long time mainly the zero-bias anomaly as well as other 
equilibrium properties, which can be extracted from the 
exact analytical solution via the Bethe ansatz approach, 
have been in the focus of the theoretical research.^^ It 
was only with the advent of nanotechnology that the in- 
vestigation of nonequilibrium properties received a boost, 
as it became possible not only to directly manufacture 
structures which are adequately described by the AIM, 
but also to investigate their nonequilibrium properties 
under well controlled parameters iLSui 

However, even in the time-independent steady-state 
case the analysis of the nonequilibrium situation turns 
out to be rather difficult. Despite a large number 
of works employing various perturbative and renor- 
malization group techniques (e.g. Refs. [10lllllll2lll3lll4 
ITEIITrilflTj ). or even attempts at solving the problem 
analytically^ there is no solution which unifies all known 
details. 

Even more difficult is the case of the "preparative" 
nonequilibrium, i.e. the time evolution of the system from 
some initial preparation towards its steady state under a 
finite external voltage bias. For the first time this prob- 
lem has been discussed in Ref. [19(, where a solution for 
the wide flat band (WFB) limit was derived. However, 
the assumption of an infinitely wide band leads to the 
rather unphysical prediction of a displacement current 
which instantaneously jumps to a finite value immedi- 
ately after switching on the tunneling. 

The transient nonequilibrium dynamics of a strongly 
interacting quantum dot which is suddenly brought into 
the Kondo regime, has been investigated using approxi- 



mative schemes . 20 ! 21 Moreover, the band structure effects 
on the time evolution of noninteracting nanoscale devices 
have been investigated in [22j. 

However, the combined effect of interaction and finite 
bandwidth, both of which can be described within the 
framework of the AIM, have not yet been considered. In 
this article we attempt to address this issue by means of 
perturbation theory in the case of weak interactions and 
Monte Carlo (MC) simulations for moderate to strong 
interactions. Such results are not only interesting for 
future experiments. In view of recent attempts to use 
the integrability methods to understand the nonequilib- 
rium properties of quantum impurity models it is im- 
portant to develop and test numerical schemes which 
are able to generate reliable results for any parameter 
constellation . 23 ' 24 ^ 

The structure of this paper is as follows. After intro- 
ducing the system under consideration in Section HU we 
start with a resonant level model which maps onto the 
noninteracting AIM. Because the corresponding Hamil- 
ton operator is quadratic in the fermionic fields, the dy- 
namics of the system can be investigated by analytic 
means at any parameter constellation even for a model 
with arbitrary band structure of the electrodes (cutoff 
schemes). The basis of our solution is the integral equa- 
tion for the impurity retarded Green's function (GF), 
which we derive next. It is then used in Section IlIII for 
the calculation of the time-dependent impurity popula- 
tion function n(t) as well as for the expectation value of 
the transient current. Here we not only consider the sim- 
plest case of a WFB structure of the leads, but also more 
realistic models taking into account bandwidth effects. 
Section IIVI is devoted to the analysis of the transient 
dynamics of an interacting system. Using perturbation 
theory in interaction strength U we identify the leading 
order effects in that limit. A treatment of arbitrary in- 
teraction strengths is best accomplished with the help of 
a dedicated nonequilibrium Monte Carlo (MC) scheme, 



2 



which is presented in Section fVl 



II. MODEL AND OBSERVABLES 

The AIM Hamiltonian usually consists of four 
contributions, 1 



In general, it is quite difficult to analyze the properties 
of the interacting Anderson model at U ^ 0, but not 
impossible. In fact, at least in equilibrium an exact ana- 
lytic solution can be derived via the Bethe ansatz^ In 
the genuine noncquilibrium, when a finite bias voltage is 
applied across the dot, the picture is far from complete 
since as yet no exact solution exists. 



H = H. 



dot 



(1) 



ifdot describes two spin-degenerate fermionic levels with 
energy A (which we later shall also call "dot"), 



H, 



dot 



= A 4d a 



(2) 



It is coupled to two fermionic continua - electrodes on the 
left and right sides. Each of these is modeled by a field 
operator ip a {x) (where a = L,R) and the corresponding 
Hamiltonians Ho[if} a ], whose precise shape we shall dis- 
cuss in a moment. The operator Ht is responsible for 
the particle exchange between the dot and the electrodes 
and is given by a simple local tunneling term of the form 



H T = E 7a tyLfr = °) d - + h - C -] 



(3) 



a=R,L a 



Finally, Hjj accounts for the interaction in the system 
and is formally implemented as an additional energy cost 
for the double occupancy of the dot level, 



Hjj = Ud\d^ d\d i 



(4) 



On the other hand, since the tunneling part Ht, being 
only quadratic in the fermionic operators, is diagonaliz- 
ablc by elementary methods and the Green's functions 
(GFs) are accessible in all parameter regimes, the non- 
interacting system [aka resonant level (RL) model] is ex- 
actly solvable by elementary means^ In this case, the 
Hamiltonian as well as expectation values separate into 
spin-up and spin-down contributions, so that throughout 
this and the next Sections we shall work with spinlcss 
operators, recovering the necessary prefactors after the 
calculations. 



In the case of the initially uncoupled dot GFs we have 
to deal with two different situations: (i) the dot level 
is empty, uq = 0, and (ii) the dot is populated by one 
electron, no = 1. Due to the simple structure of -ffdot, 
the time evolution is trivial, d(t) = d(0) exp(— iAt), im- 
mediately leading to the following matrix (Keldysh) GF 



D„(t) 



Doit) D<{t) 
D>{t) D (t) 



-iAt 



-t[e(t)(l-no)-e(-t)no] 



in 

-t[9(-t)(l-no) 



0(i)no] 



(5) 



r 



where Dq(J,) and Do(t) denote the time-ordered and anti- 
time-ordered GFs, respectively. For the retarded and 
advanced components we obtain, 



D*(t) = D (t)-D<(t) = -iQ(t)e- iAt 
D$(t) = D<(t) - D (t) = iQ(-t)e-* At 



(0) 



The Hamiltonian for the electrode electrons can generally 
be written as (for a — R, L) 



k 



(7) 



implying a trivial time evolution of the field operators. 
Due to the local tunneling assumption made in Eq. ([3]), 
coupling to the leads only involves the operator 



lp a (x = 0) = ^ ^"k ■ 



(8) 



Therefore, we only need local GFs of the band degrees of 
freedom in all subsequent calculations and can suppress 
the coordinate variable. For the retarded GF we thus 
have 



= -i0(t) / dtupiu;) 



(9) 



where we have introduced the energy-dependent density 
of states (DoS) p{ui). In a similar way, one obtains the 
full Keldysh matrix 



la(w) = i2irp(uj) 



f a ~ 1/2 f a 

-{I- fee) /a -1/2 



(10) 



Here, f a denotes the Fermi distribution function in 
the respective electrode a = L, The retarded and 
advanced components are easily retrieved, = 
—iTrp(uj) and g^{^) — [9a i 1 ^)]* — wrp(w). We would 
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like to point out that the actual dimensionality of the 
electrode disappears from the problem during the tran- 
sition from Eq. ((TJ) to ([9]) , since it is completely encoded 
in the DoS. 

The GFs of the coupled system can be found analyt- 
ically for arbitrary time dependence 7 a (£) ^ of the 
tunneling constant. From now on, we shall concentrate 
on the case of sudden switching 7 Q (t) = 7 Q O(i), where 
O(i) is the Heaviside step function. A generalization to 
arbitrary time dependence is relatively straightforward. 
The way to obtain the necessary Dyson equation is pre- 
cisely the same as in the stationary case. The result can 
be summarized in the matrix equatio n 19 i 26 (for t, t' > 0), 

/>oo />OC 

D(i,i') = T> (t-t')+ / dti \ dt 2 
Jo Jo 

x D (t-ti)£(ti-t 2 )D(t2,t , ) ) (11) 
where the generalized self-energy is defined by 

S(*)=7Lgi(*)+7«&i(*)- (12) 

The seemingly complicated structure of Eq. (TTTj) simpli- 
fies considerably for the retarded GF, 



£>«(*,*') = D R (t-t') 



(13) 



dt 2 K(t,t 2 )D H (t 2 ,t') 



where 



K(t,t 2 ) 



dt 1 D R (t-t 1 )^ R (t 1 -t 2 ) (14) 



is the kernel of the integral equation. 

The simplest physical quantity to calculate is the time- 
dependent dot population n(t) = (d* (t) d(t)) . It is con- 
venient to rewrite it in terms of the off-diagonal Keldysh 
GF, 



i(t) = -iD<(t,t) . 



(15) 



The necessary relation between this function and the al- 
ready known retarded GF is provided b y 28 ' 29 



D< = (1 + G S ) D< (1 + S G ) + G H S< G , (16) 

where products denote integration over time. This rela- 
tion is especially useful for the case of an initially empty 
dot since then Dq = and only the last term contributes. 
(A similar relation can be derived for the counterpart 
Dq , which would be useful for the initially populated 
dot.) 

Another important observable is the current through 
the system. The canonical way to calculate it is to start 
from the total particle number operator^ e.g. in the 
left electrode Ql, and to use the Heisenberg equation to 
obtain its time derivative, which is proportional to the 
current, 



Il = - 



dt 



= i[Q L ,H] 



«7i 



(17) 



The evaluation of its expectation value then leads to 
mixed correlation functions of dot and lead operators, 

hit) = d(t)) - i lL (d\t) i> L {t)) . (18) 

The expectation values entering this formula can be 
rewritten in terms of GFs. After placing the time t 
on the forward Keldysh branch and performing the con- 
tour disentanglement we obtain (we use the definition 
D K = D<+D>) 



dt 1 {[gt(t,t 1 )+g^(t,t 1 )]D A (t 1 ,t) 



+ ^(i,t 1 )[^(i 1) t)+£> ii (i 1 ,t)] 
-- [D A {t,h)+D K {t,t 1 )]gt(t 1 ,t) 



(19) 



Some products of advanced and retarded GFs vanish as 
their factors have time arguments of opposite signs. This 
simplifies the result considerably, 



i L (t) = i' L (t)+i'm, 



(20) 



where 



-y 2 r°° 

I'dt) = ^ j 

-g^{t,t l )D K {t l ,t)] 

= - 7 2 Re / dt l9 S(t,t 1 )D K (t 1 ,t), (21) 
Jo 

/>oo 

Il(t) = 7 | Re / dti D R (t, h) gf (h,t) , (22) 



and 



after using the antihermiticity of the g K and D K GFs. 
For the evaluation of these expressions it is convenient to 
use the relation D K = 2D< + D R - D A . 



III. THE NONINTERACTING CASE 

A. Wide flat band limit 

In the stationary situation, the time-translational sym- 
metry of all quantities entering Eqs. (fT3)l and (fi"4"f is 
restored and the integral equation is solved by a mere 
Fourier transformation.— In the dynamic case, the situa- 
tion is more complex. We begin with the already known 
results obtained in the approximation p(tu) = po, when 
the conduction band in the electrodes is assumed to be 
of zero curvature over an infinite range of energies. In 
this case (we concentrate henceforth on the symmetric 
coupling case 7 = 7_r, — 7l) 



9a(t) 

K(t,t 2 ) 



iirp 5(t) , 
-iirp Q 6(t) , 

_ T e -iA(t-t 2 ) e( j. . 



t2)9(i a ), 



(23) 
(24) 



where T — 2n po j 2 . As has been realized in Refs. 30lfl9l |. 
the integral equation for the retarded GF can then be 
solved by iterations. The result has the very appealing 
form 



D R (t-t') = -iQ(t-t')e 
Gathering all terms we obtain 

r 



-iA(t-i') e -r(t-t') 



(25) 



i(t) = 



2- 
1 + e 



duj \f R (u) + f L (w)] . 
-2rt_ 2e -rt cos [( w _ A)i] 



r 2 + (w - A) 2 
n sta t(l + e- 2rt ) 

' ' ' f du [/rH + /lH] 



(26) 



cos [(u - A)t] 

t 2 + {uj- A) 2 



The asymptotic, steady-state value for the population is 
given by 



r r f R (u) + f L (f>) 

"stat = — I duj 



2n J r 2 + (w - A) 2 

ii v-^ fi r 

- + — Im > * - 

2 2f ^ I 9 



(27) 



p=± 



2nT 



+ i 



pV/2 - A 
2nT 



where ^(x) denotes the psi (digamma) function. In the 
simpler case of zero temperature, it simplifies to 



"-stat 



i + -!- V arctan(pV/2 

l Z7T z — I 



A) 



(28) 



Already in Eq. (f26|) one easily identifies the tunneling 
rate T as the energy scale which governs the approach 
to the steady state. Indeed, at almost all values of other 
parameters the steady state (in which we can still have 
a finite transport current) is established after a time of 
the order r . The current which is flowing during this 
time onto and from the dot (depending on the initial con- 
dition) is to a large extent the displacement curren t 31 ! 32 
which is given by the time derivative of n(t), 



r e - rt 

Te- r * 



dn(t) 
dt 



(29) 



duj [jxM + ZiiM] 

r cos [(w - A)t] — (w — A) sin [(u - A)t] 
" r 2 + (a; - A) 2 ' 

A surprising effect is found if the dot energy is higher 
than the Fermi edges in the electrodes: on the interme- 
diate time scale A -1 , the dot population shoots over its 
asymptotic steady-state value and reaches a local maxi- 
mum despite the absence of any kind of interactions. The 
subsequent relaxation to n sta t may then be either smooth, 
or accompanied by a number of oscillations, as shown in 
Fig. [2 These are remnants of the oscillatory behavior of 
the integrand in Eq. (j2"6")l and have a period ~ A -1 . In 
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FIG. 1: Population of the noninteracting dot as a function 
of time measured in units of T" 1 for V = 0, A/F = 8 and 
different bandwidths: for e c /F = 8, 20, oo. Curves are the 
analytic results, symbols represent the MC simulation data 
(see Section fV)) . 



the opposite limit of a dot lying far below the chemical 
potentials in the leads, one does not observe related ef- 
fects. However, as the system is particle-hole symmetric, 
the analogous population oscillations are recovered for 
the initially populated dot. 

Next we investigate the time-dependent current. The 
first contribution is essentially given by the time- 
dependent n(t), 



I' L (t)= 2 (*) [l~2n(*)] . 



(30) 



du> 



leading to 

I L (t)=I L , s ^-Te~ rt J ^ {uj _ A)2+T2 

x {re- rt [f R (uj) + / L H] - rcos[(^ - A)t] [2f R {u) + 1] 
- ( W -A)sin[(^-A)t][2/ L H-1]} , 



(31) 



where the asymptotic steady-state value of the current is 

(32) 



/ r 2 / — ll^l ~ M^j 

1 ' 2tt (w- A) 2 +T 2 

I 

2 



= TGq Im > p * 



E 

P =± 



r pV/2 - A 

+ 1 



2ttT 



2-kT 



with Go = e 2 /h being the conductance quantum. The 
current I R (t) through the right contact is by symme- 
try found from I R (V, t) = — V,t). The difference 
in currents through the individual contacts describes the 
charge redistribution in the system, i.e. it must be equal 
to the displacement current 



W*) = h(V,t) + I L (-V,t) , (33) 
which is shown by inspection. The total current through 
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the AIM then becomes^, 

i(v,t) = [i L (y,t) + i a (v,t)]/2 

= [I L (V,t)-I L (-V,t)}/2. (34) 

The last two identities imply a very convenient represen- 
tation of the currents through the contacts 

7i,ii(t) = /(t)±/di S p(t)/2. (35) 

Note the instantaneous current value at t = 0, 

o f du 1 r , 

^^j2? ( M -A)HP = 2' (36) 

which is independent of both voltage and A. In fact, 
1^(0) corresponds to the current through the resonant 
level system in the case of infinitely high applied voltage. 
This unphysical, instantaneous current onset of the left 
current comes as no surprise since initially every electron 
in the band has the same probability of populating the 
empty dot, while electrons with arbitrarily high energies 
allow correspondingly fast processes. Obviously, such a 
situation can never occur in any real system due to a 
strictly finite width of the electrodes' conductance bands. 
To comply with this restraint, we proceed with analyzing 
a more realistic model with finite bandwidths. 



B. AIM with soft cutoff 

Depending on the material, the method of coupling to 
the voltage sources, and measurement apparatus the elec- 
trode band structure can be more complicated than for 
the WFB. The simplest model is a rigid flat band with a 
constant DoS po between e' c and e c , which represent the 
band bottom and upper boundary, respectively, and zero 
otherwise. A more sophisticated scheme, which is more 
physical involves soft cutoffs at e c and e' c . Mathemati- 
cally, they can be realized as 

p(u) = % ( 1 - -. ^) , (37) 

where rj is a softening parameter. From the physical point 
of view one obvious 'good' value for it would be n = T, 
which is the value we are using for numerical plots4i 
The choice of two different values, e c and e' c , is delib- 
erate. While in the equilibrium, when no bias voltage 
is applied to the system, e' c = —e c covers the relevant 
physics whatever the band filling (as long as the chemi- 
cal potentials are not too close to the band boundaries), 
the situation is more delicate when out of equilibrium. 
Usually, the finite voltage is realized by different chemi- 
cal potentials hl,r of the two electrodes. Changing them 
around the equilibrium value (we always assume the band 
to be half-filled, so that in equilibrium fiL.R — m our 
choice of zero point of energy) without shifting the band 
boundaries would imply charging of the electrodes. In 



order to avoid this, one has to shift the complete band 
along with the changed chemical potential to ensure the 
electroneutrality, e c — > e c + fiL,R, and e' c — > — e c + ^l,r- 
We shall see later that for voltages small compared to 
the band width 2e c the changes in observables vanish on 
a timescale of the order of e^ 1 . Nevertheless, in order to 
be consistent we shall keep two different values e c and e' c . 

A stronger DoS energy dependence is expected for a 
system strongly coupled to its environment. In some 
cases even a DoS which vanishes at the Fermi level 
may emerge. Two notable situations are a system in 
the Coulomb blockade regime or a strongly interacting 
low-dimensional conductor in the Luttinger liquid phase. 
These are, however, systems with effectively interacting 
electrodes whose treatment we postpone to a follow-up 
publication. From now on, we would like to concentrate 
on the DoS (J37J). 

In order to solve the Dyson equation (fT"3"]) it is more 
convenient to have the retarded band GF in the time 
representation. According to the prescription ^ we find 

g— ie' c t _ g—itct 

™*W [e «^-l]sinh(^) ' (38) 

which is regular towards the limit t = 0. From this rela- 
tion one can easily read off the retarded self-energy using 
Eq. (fT2)) . The corresponding integral equation kernel (jT4j) 
turns out to depend only on the time difference (t — t2), 

rTe" lA(t "* 2) 
KCt-h) = -i ~, ttt (39) 

x dr e tAr . 

Jo sinh(7r?7r) 

In fact, the last integral can be expressed in terms of hy- 
pergeometric functions. However, from a numerical point 
of view, it is more convenient to work with the integral 
(|39|) directly. In fact, writing down the equations for the 
retarded GFs in the steady state case and in the case of 
the sudden switching of tunneling one immediately real- 
izes that they are identical in the relevant time domain. 
In order to calculate the time-dependent population func- 
tion, one still can use the formula p^|) . The necessary 
off-diagonal self-energy is given by 

VTp-^/ T ^ 

s<(*) = - r. u rp.\ E ( 4 °) 

2sinh(7rTt) 

e -ip, z t-^/T 

X [(e-w/T _ e - e i/T)( e - Mi /r _ e -« c /T) 

e -ie' c t-e'JT 

~ (e-</ T - e-«WT)( e -M</T _ e -</T) 

e -ie c t-e c /T 

+ {e-</ T - e-WT)( e -WT _ e -<W)_ ' 

With the prerequisites (|3"9")l and (|4"U)) the calculation of 
the time evolution n(t) as well as of the currents is a 



FIG. 2: DQMC data for the left, right, and total current 
lL,n(t), and I(t), respectively, for U = A = 0, V = 20r, and 
T — OAT. Circles, triangles down, and triangles up refer to 
cutoff energies of e c = 20F, 40r, and 100r, respectively. The 
inset shows the difference in I(t) due to charge neutrality. 
The solid lines refer to the current calculated in the WFB 
limit via Eq. ((3Tj) . 



FIG. 3: DQMC data for the displacement current idisp(t) 
(open symbols) and results according to the perturbation the- 
ory (lines) for V = 0, e c = 20F, and T = OAT. Circles, 
squares, diamonds, and triangles (solid, dashed, dotted, and 
dotted-dashed lines) refer to U = -2A = 0, T, 2V, and 4F, 
respectively. 



IV. PERTURBATION THEORY IN 
INTERACTION 



standard numerical task. The results of the calculations 
are presented in Figs. Q] and [21 The most drastic dif- 
ferences between the WFB model and that with a finite 
bandwidth are found in the short time behavior of the 
current. In contrast to the WFB prediction the instan- 
taneous value of currents through individual contacts is 
strictly zero. Moreover, the slope (time derivative) of 
lL,R,(i) starts at zero rather than being finite. These dif- 
ferences eventually vanish after a timescale of the order 
e" 1 , so that, as expected, the correspondence between the 
two calculation schemes improves. However, the actual 
current behavior becomes more oscillatory and prevents 
reliable simulations for too large e c /T > 50. Contrary to 
the lL,n{t) currents, the total current through the con- 
striction is not only free of oscillations but also shows a 
far better agreement with the WFB model. We conclude 
that the finite bandwidth effects are contained almost 
completely in the displacement component of the current 
(|2"5)) . Its behavior is plotted in Fig. [3] 

Furthermore, we find a rather small difference between 
the results for systems which preserve and neglect the 
electroneutrality (see inset of Fig. \2§ , which only exists 
for times ~ e^ 1 and vanishes almost completely in the 
steady state. We would like to point out that the max- 
imal deviation depicted in Fig. [5] is achieved for volt- 
ages half as large as the bandwidth. It is highly unlikely 
that such a situation can ever be realized in experiments, 
where the maximal voltages very seldom exceed 5% of 
2e c . Therefore from now on we refrain from implement- 
ing the electroneutrality requirement in our analysis. 



As a next step, we investigate the change of the dot 
transient dynamics due to the finite Coulomb repulsion, 
which is described by the term Q in the Hamiltonian. 
In the regime where U is small compared to the other 
energy scales, we can employ a perturbative expansion 
which we truncate after the first order. Note that as the 
interaction involves electrons of opposite spins, we have 
to keep track of the spin indices henceforth. 

In order to calculate the time-dependent dot popula- 
tion n a (t), we start from Eq. (fT5")) and expand the dot GF 
to first order in U. Discarding all disconnected diagrams, 
this leads to 

D^(t,t') = U [ dsn*(s)D a (t,s)D a (s,t'). (41) 
Jc 

The superscript denotes the first order in U, while the 
GFs on the right hand side and the particle density n s (s) 
are the respective functions for the noninteracting case. 
The time integration runs along the Keldysh contour C. 
The lesser GF can be expressed in terms of the advanced 
and retarded GFs and reads 

/oo 
ds n^s)[D^(t,s)D<(s,t') 
-CO 

+ D<(t,s)D£(s,t')]. (42) 

Both components can be combined by using the complex 
conjugation properties 

D$(t,t) = [D?(t',t)r, 

D<{t,t>) = -[D<{t',t)}*. (43) 



Initially, before the coupling to the leads is switched 
on, the dot is assumed to be empty, which means that 
n s (s) = for s < 0. Therefore, we can rewrite the first 
order correction to the dot occupation number as 

n£>(t) = 2 UIm[r a (t)} , (44) 



r a (t)= dsn 9 (s)D^{t,s)D<(s,t) (45) 



where 



depends only on properties of the nonintcracting system 
and is thus accessible. In the WFB limit, this calculation 
can mostly be done analytically using the functions given 
in Eqs. (53) and (25]). 

In order to keep the derivation simple, we shall inves- 
tigate the case A = in equilibrium and at zero temper- 
ature (V = T = 0). The change in the asymptotic value 
can most easily be accessed, because the usual pertur- 
bation theory in the frequency domain can be employed. 
One finds the following correction to the dot occupation, 



D a {u)) - D a {u) 



(46) 

Using the well-known expressions for the dot GFs D l J (lo) , 
and the fact that the unperturbed stationary dot occu- 

1/2, one finds 



pation is given by n^ tat 



U 
2nT 



(47) 



which is a simple reduction of the stationary value due to 
the Coulomb repulsion. In a next step, we shall find out 
how this stationary value is approached. For this pur- 
pose, we evaluate the function (|4"5)) which contains the 
unperturbed time-dependent occupation number. For 
our choice of parameters, this function can be read off 
from Eq. (j26j) : 



4 0) w^(i- e - 2rt ) 



(48) 



Hence, in equilibrium the dot occupation without 
Coulomb interaction saturates exponentially on a time- 
scale r _1 . Moreover, we need the lesser dot GF which 
can be derived from Eq. (|16[) . One finds 



D<(s,t) 



IT 



II 



dtu 



r 2 



e -Vs e iut 



e -rt e -ius 



(49) 

-T{s+t) 



For A = 0, the retarded GF (|25|) becomes purely imagi- 
nary and the only imaginary part in Eq. (j45j arises from 
the lesser GF D<(s,t). Therefore, we only need to eval- 
uate the imaginary part of the w-integrals. With the 
definitional 
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FIG. 4: Time-dependent dot occupation for different interac- 
tion strengths U (V = 0, A = -(7/2, e c /Y = 10, T = 0.2r). 
Curves show the results from the first order perturbation cal- 
culation in U . The symbols represent the MC data. 



where Shi (a;) and Chi (a;) denote the hyperbolic sine and 
cosine integral functions, respectively, we obtain 



«W(t) 



2UT 



4 0) ( S )e- r C- 



(51) 



x [z(t — s) 



z(t) 



This integral is evaluated numerically for arbitrary t. As 
the integrand is regular towards s,t — > 0, one concludes 
that the time derivative vanishes for small times and the 
correction due to the Coulomb interaction only sets in 
gradually. This is understandable since we assumed the 
dot to be initially unoccupied and the Coulomb interac- 
tion can only have an effect once a finite population has 
been established on the dot. 

The previous calculations for V = A = T = clearly 
present an oversimplified picture. Although the qualita- 
tive statements remain correct, additional energy scales 
due to finite temperature, asymmetry, voltage, and band- 
width render the whole picture more complicated. In 
these cases, however, an all- numerical scheme has to be 
used. We chose to solve the integral equation (fl"3|) by 
discretizing the time axis and thus translating it into a 
matrix equation. 

In order to investigate the limit of the approximation 
in U, we compared the perturbative corrections to MC 
results which are exact for any interaction strength. In 
Fig. [H we plot the time-dependent occupation proba- 
bility for a single spin orientation in the noninteracting 
model and for a relatively small interaction U. While for 
small times, the graphs coincide within numerical accu- 
racy, deviations become visible after a time of order 
This is not surprising because in order for the interac- 
tion term to be fully operational a finite dot population 
is necessary. This process requires a time of the order 

r- 1 . 

In order to investigate the interaction effect more 



closely, we plot in Fig. [5] the interaction correction for 
several values of U. One notices a good agreement up to 
times of order T _1 even for interaction strengths as large 
as U = r, which is actually beyond the range of the per- 
turbation theory. In this regime, even an expansion up 
to second order in U, albeit feasible in principle, would 
not lead to a more reliable result. For this reason, we 
only consider the first order in U. 

While the effect of the Coulomb interaction on the dis- 
placement current can be calculated by differentiating 
Eq. (fSTj) with respect to time, we still have to investigate 
the time dependence of the average current. As we ar- 
gued previously, its dependence on the electronic band- 
width is very small, so we shall do this analysis in the 
WFB limit. The calculation starts again from Eqs. (j2"Tj) 
and (|22p . where we have to use the first order expansions 
of the dot GFs. Hence the expression for the current 
across, say the left lead, is given by a sum of two terms 
1^ = I'l (t) + l'^ X \t), where the first one is propor- 
tional to the dot occupation 

l'£{t) = TQ[t)n£Xt) . 



(52) 



The derivation of the second contribution involves the 
complete set of dot GFs, D K , D R and D A . A straightfor- 
ward calculation then leads to the following result which 
is valid at zero temperature, 



r"(l) 
l Lcr 



(*) = 



eUT 



ds 



,-T(t-s) 



(53) 



X {- - cos[|V72 - A\(t - «)]} J* ds' nf\s'). 

This integral can be calculated numerically using the 
known zero-order dot population n§\t). The result is 
compared with the MC data in Fig. [6l In contrast to the 
dot occupation, the agreement between perturbation the- 
ory and numerically exact MC simulations for the current 
degrades rapidly towards stronger interaction, see Fig. [6] 
for U = 8r. The steady state value of the current is con- 
siderably smaller than that for weak interactions. How- 
ever, this does not contradict the conventional Kondo 
picture, where the conductance is enhanced, since in our 
case the voltage is higher than the interaction strength. 
Another feature is the extremely weak oscillations in re- 
lation to the small U case, which can be interpreted as a 
precursor of the Kondo physics as discussed in [2l[ . 



V. MONTE CARLO SIMULATIONS AT 
ARBITRARY INTERACTION STRENGTH 

To go beyond the first order perturbation theory, we 
employ a numerically exact diagrammatic Monte Carlo 
scheme. The method is a generalization of the algorithm 
proposed in Ref. (35| . or - equivalently - the Keldysh 
implementation of the diagrammatic impurity solver of 
Ref. [36j. Here, we only outline the basic concepts of the 
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FIG. 5: Change of the time-dependent dot occupation due 
to interactions, 5n a (U) = n a (U) — n a {TJ = 0), for several 
values of the interaction strength. V = 0, A = —(7/2, e c /Y — 
10, T = 0.2r. Curves show the results from the first order 
perturbation calculation in U. The symbols represent the MC 
data. 



algorithm and refer to Refs. 35,36] for further details; 
the specific aspects of dealing with the electron-electron 
correlations in the real-time simulations will be presented 
in a forthcoming publication. 

The main idea is to evaluate the expectation value 
(0(t)) = Tr dotilead [^ ot ' lead e itH Oe- rtff ] by stochastically 
sampling a perturbation expansion in the tunneling term 
Ht- To this end, we employ an interaction represen- 
tation in which the time evolution along the Kadanoff- 
Baym contour — > t — > is determined by _ffi oc + Ho = 
Hdot + H\j + Hq and rewrite the time evolution operators 



,±itH 



as (anti-)time ordered exponentials 



(0(*)> 



Tr, 



dot. lead 



Te -ifid,H T (.s)] (54) 



With 0(s) = e it(,H loc +H ) Oe -it(H lac +H ) ( and Ht ( s } ac _ 

cordingly). The exponentials are then expanded into a 
power series, which allows to trace out the electrode de- 
grees of freedom in an exact manner. At perturbation 
order N (for given spin), this yields a determinant of 
an N x N matrix whose elements are determined by the 
self-energy and the times at which the tunneling events 
occur. The configuration space consists of all possible se- 
quences of dot creation and annihilation operators on the 
Kadanoff-Baym contour, and the MC sampling proceeds 
through local updates of these operator sequences (inser- 
tion/removal of pairs of creation and annihilation oper- 
ators, or shifts of the operator positions). We use both 
a continuous-time implementation (CTQMC) as well as 
one which utilizes a discretization scheme for the real- 
time axis to speed up the sampling process (DQMC). 
While the continuous-time approach is completely free 
of systematic errors, we kept discretization effects be- 
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FIG. 6: DQMC data for the total current I(t) (symbols) and 
results according to Eqs. (|52[l and (|53[l (lines) for V = 10r, 
T = 0, and e c = I0F and 40r (open and filled symbols, respec- 
tively). A WFB was assumed for the first order perturbation 
calculation. Circles, squares, diamonds, and triangles (facing 
down and up) refer to A = -17/2 and U = 0, F, 2T, 4T, and 
8F, respectively, while solid, dashed, and dotted lines refer to 
U = 0, T, and 2r. 



low statistical errors in the discrete time implementation 
as well, by using a fine mesh of 10 3 points along the 
real-time axis. From the simulation, we obtain the time- 
dependent dot population n(t) as well as Ih,ii(t), Idi sp (t), 
and the total current I(t). 

While this simulation approach can handle arbitrary 
interaction strengths, it suffers from a dynamical sign 
problem™ which becomes severe at long times or for large 
bandwidth. We find that the error bars grow exponen- 
tially with average perturbation order, and thus expo- 
nentially with time t. In the noninteracting model, which 
factorizes into spin-up and -down components, only one 
spin species needs to be simulated. This reduces the av- 
erage perturbation order by a factor of two and allows 
us to simulate a time interval which is about twice as 
long as in an interacting model (in contrast to imaginary 
time, the perturbation order is essentially independent of 
interaction strength) . Our simulations of the interacting 
AIM can reach the stationary state for certain parameters 
(cf. Fig. IHJ), but not yet for in the general case. We note 
in passing that in principle, the inclusion of a phonon 
background coupling to the interacting dot is straight- 
forward within the path- integral framework of Ref . [35| , 
or using the method proposed in Ref. [371 ] . 



VI. DISCUSSION AND CONCLUSIONS 

After switching on the tunneling, the dynamics of the 
AIM exhibits a surprisingly rich transient behavior. The 
reason can be found in the abundance of different energy 
scales. While in the steady state, all parameters like 
temperature T, voltage V, contact transparency T, and 
dot parameters - the bare energy A and the interaction 
strength U - are known to be decisive for the stationary 
values of the transport current as well as the population 
probability of the dot, in the opposite limit of intermedi- 
ate and especially short times the dynamics is dominated 
by the influence of T and A. In fact, in both the interact- 
ing and the noninteracting case the typical timescale at 
which the steady state is reached is of the order T -1 . In 
the noninteracting case it is clearly visible that the ap- 
proach towards steady state is exponential, see Eqs. ([^rj]) 
and (|3ip . On the other hand, a nonzero detuning A or 
a finite band cutoff lead to superimposed oscillations of 
the observables' time evolution. 

The most interesting behavior is encountered at very 
short timescales. It turns out that the rather simple 
WFB theory fails to give meaningful results here, pre- 
dicting e.g. an instantaneous finite value for the currents 
through the individual contacts (|36|) . This unphysical 
picture can only be corrected by considering a more re- 
alistic model for the electrodes featuring a finite band- 
width e c , which then dominates the short time dynamics 
of the system. It slows down the onset of the current and 
dot population and thus quenches their time derivative 
to much smaller values than for an infinite cutoff. Only 
after a timescale of the order of e" 1 do current and dot 
population approach the values of the WFB model. 

However, this comes as no surprise since for a system 
with a finite range of allowed excitations W (in our case 
the electrodes with W ~ e c ), the uncertainty principle 
demands that the reaction to any instantaneous pertur- 
bation (switching on of tunneling) has to take place on 
a finite timescale ~ W^ 1 . Furthermore, the MC simula- 
tions reveal fast oscillations in the currents through the 
individual contacts, whose wavelength and amplitude de- 
creases with increasing e c . 

In contrast to the currents through the electrodes 
and the diplacement current, however, the total current 
through the system is only weakly affected by the band- 
width. According to the numerical results, even for band- 
widths only twice as large as the voltage, the current 
follows the analytical results for the WFB with high ac- 
curacy. This stems from the fact that the displacement 
current (|29[) monitors the redistribution of charge across 
the electrodes, but it is not responsible for any net charge 
flow; on the other hand, the cutoff effects are almost com- 
pletely due to the charge redistribution. Therefore, it 
appears natural to expect the total current to be only 
weakly affected by the explicit value of e c . We find virtu- 
ally no influence on the long-time dynamics for realistic 
voltage/bandwidth quotients. The maximal deviations 
with respect to the WBF limit is again achieved during 
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times ~ e^ 1 . 

For an initially empty dot the effects of a finite 
Coulomb interaction become visible only after a timescale 
r . This can be rationalized by observing that the same 
timescale is necessary to build up a dot population large 
enough to be affected by electronic correlations. Further- 
more, we find that the quality of the approximation by 
lowest order perturbation expansion is remarkably good 
(even for intermediate U) up to tT rj 1, which is about 
a factor of three smaller than the time required to reach 
the steady state. 

This shows that similarly to the Kondo case of 
Ref. [38|, the full interaction effects take a time of sev- 
eral F _1 to develop. Another interaction effect is the 
suppression of both the dot population and the current. 
While the first effect is quite natural, the second one 
is seemingly at odds with the common wisdom that at 
sufficiently low temperatures, due to the Kondo effect, 
the transport properties of the system must approach 
the unitary limit of a perfectly resonant level, see e.g. 
Ref. [39[ . Our results do not allow to see this kind of 
physics for two reasons: (i) the steady state is not yet 
fully established, (ii) the applied voltage is rather large 
and therefore has the potential to destroy the Kondo ef- 



fect even in the steady state regime. 

To conclude, we presented a theoretical treatment of 
transient effects in an AIM biased with a finite voltage 
after a sudden switching on of the tunneling. Using ex- 
act analytical solutions and perturbation theory as well 
as dedicated numerical schemes (MC) we identified dif- 
ferent regimes in the time evolution of the currents and 
the dot's population probability and related them to the 
parameters of the system. Special attention has been 
paid to the influence of electron-electron interactions on 
the dot and the bandwidth of the electrodes. 
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We assume the electrodes to be large enough in order to 
have a meaningful concept of temperature. 
At finite temperature the single particle levels can be con- 
sidered as broadened, having the typical width of T. That 



is why it is natural to assume the band boundaries to be 
widened with rj ~ T. 

In imaginary time, the diagrammatic Monte Carlo simula- 
tion of the AIM does not encounter a sign problem. All the 
(fermionic) sign cancelations between crossing and non- 
crossing diagrams can be absorbed into the determinant 
of self-energies (hybridization functions), and the imagi- 
nary time evolution e~ HlacT gives a real contribution to 
the weight. In real-time diagrammatic Monte Carlo, both 
the self-energies and the time evolution operator are com- 
plex. 



